# SimesJahnSI.R
# Supporting Information
# "The Consequences of Medicaid Expansion Under the Affordable Care Act for Police Arrests"
# Last updated by JS: 7/13/21

library(dplyr)
library(ggplot2)
library(sandwich)
library(cowplot)
library(MatchIt)
library(optmatch)
library(tidyr)

##-----------------------------------------------------------------------
## Reading in main data, Sandwich SEs function
##-----------------------------------------------------------------------

load("SimesJahn.RData")

ses<- function(model) { 
  cov <- vcovHC(model, type = "HC0")
  std.err <- sqrt(diag(cov))
  q.val <- qnorm(0.975)
  r.est <- as.data.frame(cbind(
    Estimate = coef(model)
    , "Robust SE" = std.err
    , z = (coef(model)/std.err)
    , "Pr(>|z|) "= 2 * pnorm(abs(coef(model)/std.err), lower.tail = FALSE)
    , LL = coef(model) - q.val  * std.err
    , UL = coef(model) + q.val  * std.err
  ))
  return(r.est)
}

##---------------------------------------------------------------------------
## Setting up Propensity Score Matching Procedure for SI
##---------------------------------------------------------------------------

# reshaping in wide file for matching
w <-analysis %>% 
              gather(variable,value,-(c(FIPS:ctyname,urbanicity:studycounty,expansion:nonexp))) %>% 
              unite(temp, Year, variable) %>% 
              spread(temp, value)
w <-w[!(is.na(w$"2011_all")==TRUE),]
colnames(w)<-c("FIPS","state","ctyname","urbanicity","urban","studycounty","expansion","treatment","Expansion14",
                           "Expansion15","Expansion16","nonexp","all2011","arr2011","black2011","ci2011","countypop2011","coverage_indicator2011",
                           "crimer2011","death2011","drr2011","drug2011","exeduc2011","exsocial2011","ldeath2011","lowlevel2011","medage2011",
                           "poverty2011","unemp2011","violent2011","all2012","arr2012","black2012","ci2012","countypop2012","coverage_indicator2012",
                           "crimer2012","death2012","drr2012","drug2012","exeduc2012","exsocial2012","ldeath2012","lowlevel2012","medage2012",            
                           "poverty2012","unemp2012","violent2012","all2013","arr2013","black2013","ci2013","countypop2013","coverage_indicator2013",
                           "crimer2013","death2013","drr2013","drug2013","exeduc2013","exsocial2013","ldeath2013","lowlevel2013","medage2013",
                           "poverty2013","unemp2013","violent2013","all2014","arr2014","black2014","ci2014","countypop2014","coverage_indicator2014",
                           "crimer2014","death2014","drr2014","drug2014","exeduc2014","exsocial2014","ldeath2014","lowlevel2014","medage2014","poverty2014",            
                           "unemp2014","violent2014","all2015","arr2015","black2015","ci2015","countypop2015","coverage_indicator2015","crimer2015",
                           "death2015","drr2015","drug2015","exeduc2015","exsocial2015","ldeath2015","lowlevel2015","medage2015","poverty2015",          
                          "unemp2015","violent2015","all2016","arr2016","black2016","ci2016","countypop2016","coverage_indicator2016","crimer2016",
                           "death2016","drr2016","drug2016","exeduc2016","exsocial2016","ldeath2016","lowlevel2016","medage2016","poverty2016","unemp2016","violent2016" )

pre_cov <- c('medage2011', 'medage2012', 'medage2013', 
             'black2011', 'black2012', 'black2013', 
             'poverty2011','poverty2012','poverty2013',
             'unemp2011','unemp2012','unemp2013',
             'exsocial2011', 'exsocial2012','exsocial2013', 
             'exeduc2011','exeduc2012','exeduc2013',
             'ldeath2011','ldeath2012','ldeath2013')

# Propensity score matching
m_ps <- glm(treatment ~ medage2011 + medage2012 + medage2013
            + black2011 + black2012 + black2013 
            + poverty2011 + poverty2012 + poverty2013
            + unemp2011 + unemp2012 + unemp2013
            + exsocial2011 + exsocial2012 + exsocial2013
            + exeduc2011 + exeduc2012 + exeduc2013
            + ldeath2011 + ldeath2012 + ldeath2013,
            family = binomial(), data = w)
summary(m_ps)

prs_df <- data.frame(pr_score = predict(m_ps, type = "response"),
                     treatment = m_ps$model$treatment,
                     fips = m_ps$data$FIPS)

# Integrating PS into main data for plotting S3 Fig
idind <- match(analysis$FIPS,prs_df$fips)
analysis$ps <- ps <- prs_df$pr_score[idind]

# Matching algorithm
prenomiss <- w %>%  # MatchIt does not allow missing values
  select(FIPS, treatment, medage2011, medage2012, medage2013, black2011, black2012, black2013,
         poverty2011, poverty2012, poverty2013, unemp2011, unemp2012, unemp2013,
         exsocial2011, exsocial2012, exsocial2013, exeduc2011, exeduc2012, exeduc2013, 
         ldeath2011, ldeath2012, ldeath2013, one_of(pre_cov)) %>%
  na.omit()

pre_mod_match <- matchit(treatment ~ medage2011 + medage2012 + medage2013
                         + black2011 + black2012 + black2013 
                         + poverty2011 + poverty2012 + poverty2013
                         + unemp2011 + unemp2012 + unemp2013
                         + exsocial2011 + exsocial2012 + exsocial2013
                         + exeduc2011 + exeduc2012 + exeduc2013
                         + ldeath2011 + ldeath2012 + ldeath2013,
                         data=prenomiss,
                         method = "nearest", replace=FALSE)

dta_m <- match.data(pre_mod_match)

# Matching
idind <- match(analysis$FIPS,dta_m$FIPS)
analysis$distance <- distance <- dta_m$distance[idind]
analysis$pssample <- ifelse(is.na(analysis$distance), 0,1)
matched_sample <-data.frame(analysis[analysis$pssample==1,])

##-----------------------------------------------------------------------
## S2 Coverage Indicator Table
##-----------------------------------------------------------------------

#S2 Table
analysis %>% group_by(ci) %>% summarise(arrests=mean(arr,na.rm=TRUE),nexp=mean(nonexp,na.rm=TRUE),exp14=mean(Expansion14,na.rm=TRUE), 
                                      exp15=mean(Expansion15,na.rm=TRUE),exp16=mean(Expansion16,na.rm=TRUE),black=mean(black,na.rm=TRUE),
                                      mage=mean(medage,na.rm=TRUE),pov=mean(poverty,na.rm=TRUE),unemp=mean(unemp,na.rm=TRUE),
                                      urb=mean(urban,na.rm=TRUE),welf=mean(exsocial,na.rm=TRUE),edu=mean(exeduc,na.rm=TRUE),
                                      pop=mean(countypop,na.rm=TRUE), opioidd=mean(death, na.rm=TRUE))

summary(aov(arr ~ ci, data=analysis))
summary(aov(countypop ~ ci, data=analysis))
summary(aov(black ~ ci, data=analysis))
summary(aov(medage ~ ci, data=analysis))
summary(aov(poverty ~ ci, data=analysis))
summary(aov(unemp ~ ci, data=analysis))
summary(aov(exsocial ~ ci, data=analysis))
summary(aov(exeduc ~ ci, data=analysis))
summary(aov(death ~ ci, data=analysis))

t1<-table(analysis$urban, analysis$ci)
chisq.test(t1)

t1<-table(analysis$nonexp, analysis$ci)
chisq.test(t1)
t1<-table(analysis$Expansion14, analysis$ci)
chisq.test(t1)
t1<-table(analysis$Expansion15, analysis$ci)
chisq.test(t1)
t1<-table(analysis$Expansion16, analysis$ci)
chisq.test(t1)

##-------------------------------------------------------------------------
## Placebo control models
##-------------------------------------------------------------------------
library(MASS)

load("SimesJahn_placebo.RData")

result_all <- as.data.frame(matrix(ncol = 7, nrow = 11)) 
colnames(result_all)<-c("year","RR","SE","Z","P","LL","UL")
for (i in 1:11) {
  analysisp$placebo<-ifelse(analysisp$Year>=2000+i,1,0)
  
  mp1<-glm.nb(all ~ as.factor(State) + as.factor(Year)+ as.factor(placebo)*nonexp
              + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
              + offset(log(countypop)), data=analysisp)
  result_all[i,1]<- 2000+i
  result_all[i,2]<- ses(mp1)[68,1]
  result_all[i,3]<- ses(mp1)[68,2]
  result_all[i,4]<- ses(mp1)[68,3]
  result_all[i,5]<- ses(mp1)[68,4]
  result_all[i,6]<- ses(mp1)[68,5]
  result_all[i,7]<- ses(mp1)[68,6]
}
result_all$model<-"All"

result_vio <- as.data.frame(matrix(ncol = 7, nrow = 11)) 
colnames(result_vio)<-c("year","RR","SE","Z","P","LL","UL")
for (i in 1:11) {
  analysisp$placebo<-ifelse(analysisp$Year>=2000+i,1,0)
  mp2<-glm.nb(violent ~ as.factor(State) + as.factor(Year)+ as.factor(placebo)*nonexp
              + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
              + offset(log(countypop)), data=analysisp)
  result_vio[i,1]<- 2000+i
  result_vio[i,2]<- ses(mp2)[68,1]
  result_vio[i,3]<- ses(mp2)[68,2]
  result_vio[i,4]<- ses(mp2)[68,3]
  result_vio[i,5]<- ses(mp2)[68,4]
  result_vio[i,6]<- ses(mp2)[68,5]
  result_vio[i,7]<- ses(mp2)[68,6]
}
result_vio$model<-"Violent"

result_drug <- as.data.frame(matrix(ncol = 7, nrow = 11)) 
colnames(result_drug)<-c("year","RR","SE","Z","P","LL","UL")
for (i in 1:11) {
  analysisp$placebo<-ifelse(analysisp$Year>=2000+i,1,0)
  mp3<-glm.nb(drug ~ as.factor(State) + as.factor(Year)+ as.factor(placebo)*nonexp
              + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
              + offset(log(countypop)), data=analysisp)
  result_drug[i,1]<- 2000+i
  result_drug[i,2]<- ses(mp3)[68,1]
  result_drug[i,3]<- ses(mp3)[68,2]
  result_drug[i,4]<- ses(mp3)[68,3]
  result_drug[i,5]<- ses(mp3)[68,4]
  result_drug[i,6]<- ses(mp3)[68,5]
  result_drug[i,7]<- ses(mp3)[68,6]
}
result_drug$model<-"Drug"

result_low<- as.data.frame(matrix(ncol = 7, nrow = 11)) 
colnames(result_low)<-c("year","RR","SE","Z","P","LL","UL")
for (i in 1:11) {
  analysisp$placebo<-ifelse(analysisp$Year>=2000+i,1,0)
  mp5<-glm.nb(lowlevel ~ as.factor(State) + as.factor(Year)+ as.factor(placebo)*nonexp
              + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
              + offset(log(countypop)), data=analysisp)
  result_low[i,1]<- 2000+i
  result_low[i,2]<- ses(mp5)[68,1]
  result_low[i,3]<- ses(mp5)[68,2]
  result_low[i,4]<- ses(mp5)[68,3]
  result_low[i,5]<- ses(mp5)[68,4]
  result_low[i,6]<- ses(mp5)[68,5]
  result_low[i,7]<- ses(mp5)[68,6]
}
result_low$model<-"Low Level"

mods<-rbind(result_all,result_vio,result_drug,result_low)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$RR<-exp(mods$RR)
mods$Change<-(mods$RR-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100
mods$year<-as.factor(mods$year)

#pdf("S1fig.pdf")
p1<-ggplot(mods[mods$model=="All",], aes(x=year, y=Change))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL))+
  coord_cartesian(ylim=c(-12,12))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("All")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p1
p2<-ggplot(mods[mods$model=="Violent",], aes(x=year, y=Change))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL))+
  coord_cartesian(ylim=c(-12,12))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("Violent")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p2
p3<-ggplot(mods[mods$model=="Drug",], aes(x=year, y=Change))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL))+
  coord_cartesian(ylim=c(-12,12))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("Drug")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p3
p4<-ggplot(mods[mods$model=="Low Level",], aes(x=year, y=Change))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL))+
  coord_cartesian(ylim=c(-12,12))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("Low Level")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p4
#dev.off()

##-------------------------------------------------------------------------
## Sensitivity Analysis of DID Models: Crime Adjusted
##-------------------------------------------------------------------------

m1c<-glm.nb(all ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath + crimer
            + offset(log(countypop)), data=analysis)
m2c<-glm.nb(violent ~ as.factor(state) + as.factor(Year)+ as.factor(expansion) 
            + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath + crimer 
            + offset(log(countypop)), data=analysis)
m3c<-glm.nb(drug ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath + crimer
            + offset(log(countypop)), data=analysis)
m4c<-glm.nb(lowlevel ~ as.factor(state) + as.factor(Year)+ as.factor(expansion) 
            + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath + crimer
            + offset(log(countypop)), data=analysis)

mod1<-ses(m1c)[52:54,]
mod2<-ses(m2c)[52:54,]
mod3<-ses(m3c)[52:54,]
mod4<-ses(m4c)[52:54,]

mods<-rbind(mod1,mod2,mod3,mod4)
mods$Estimate<-exp(mods$Estimate)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$model<-c(rep("All",3),rep("Violent",3),rep("Drug",3),rep("Low Level",3))
mods$year<-c(rep(c(1,2,3),4))
mods$Change<-(mods$Estimate-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100

write.csv(mods, file="results_crimeadj.csv")

##-------------------------------------------------------------------------
## Sensitivity Analysis of DID Models: State-Level Models
##-------------------------------------------------------------------------

stanalysis <- data.frame(read.csv("SimesJahn_state.csv", header=TRUE, stringsAsFactors =FALSE))

# Note: we exclude the log opioid mortality rate due to non-convergence

# State Level Model: All arrests
m1s<-glm.nb(all ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + black + poverty + unemp + exsocial + exeduc
            + offset(log(totpop)), data=stanalysis)
m2s<-glm.nb(violent ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + black + poverty + unemp + exsocial + exeduc 
            + offset(log(totpop)), data=stanalysis)
m3s<-glm.nb(drug ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + black + poverty + unemp + exsocial + exeduc 
            + offset(log(totpop)), data=stanalysis)
m4s<-glm.nb(lowlevel ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + black + poverty + unemp + exsocial + exeduc 
            + offset(log(totpop)), data=stanalysis)

mod1<-ses(m1s)[52:54,]
mod2<-ses(m2s)[52:54,]
mod3<-ses(m3s)[52:54,]
mod4<-ses(m4s)[52:54,]

mods<-rbind(mod1,mod2,mod3,mod4)
mods$Estimate<-exp(mods$Estimate)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$model<-c(rep("All",3),rep("Violent",3),rep("Drug",3),rep("Low Level",3))
mods$year<-c(rep(c(1,2,3),4))
mods$Change<-(mods$Estimate-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100
write.csv(mods, file="results_state.csv")

##-------------------------------------------------------------------------
## Sensitivity Analysis of DID Models: Propensity Score Matching Models
##-------------------------------------------------------------------------

# Matching Model: All arrests
m1m<-glm.nb(all ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + urbanicity + unemp
            + medage + black + poverty + exsocial + exeduc + ldeath
            + offset(log(countypop)), data=matched_sample)

# Matching Model: Violent arrests
m2m<-glm.nb(violent ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + urbanicity + unemp
            + medage + black + poverty + exsocial + exeduc + ldeath
            + offset(log(countypop)), data=matched_sample)

# Matching Model: Drug arrests
m3m<-glm.nb(drug ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + urbanicity + unemp
            + medage + black + poverty + exsocial + exeduc + ldeath
            + offset(log(countypop)), data=matched_sample)

# Matching Model: Low-level arrests
m4m<-glm.nb(lowlevel ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
            + urbanicity + unemp
            + medage + black + poverty + exsocial + exeduc + ldeath
            + offset(log(countypop)), data=matched_sample)

mod1<-ses(m1m)[51:53,]
mod2<-ses(m2m)[51:53,]
mod3<-ses(m3m)[51:53,]
mod4<-ses(m4m)[51:53,]

mods<-rbind(mod1,mod2,mod3,mod4)
mods$Estimate<-exp(mods$Estimate)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$model<-c(rep("All",3),rep("Violent",3),rep("Drug",3),rep("Low Level",3))
mods$year<-rep(c(1,2,3),4)
mods$Change<-(mods$Estimate-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100
write.csv(mods, file="results_psm.csv")

##-------------------------------------------------------------------------
## Sensitivity Analysis of DID Models: Coverage Indicator Adjusted
##-------------------------------------------------------------------------
ne0<-analysis[(analysis$coverage_indicator!=0),]

m1<-glm.nb(all ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=ne0)
m2<-glm.nb(violent ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath 
           + offset(log(countypop)), data=ne0)
m3<-glm.nb(drug ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=ne0)
m4<-glm.nb(lowlevel ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=ne0)

mod1<-ses(m1)[52:54,]
mod2<-ses(m2)[52:54,]
mod3<-ses(m3)[52:54,]
mod4<-ses(m4)[52:54,]

mods<-rbind(mod1,mod2,mod3,mod4)
mods$Estimate<-exp(mods$Estimate)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$model<-c(rep("All",3),rep("Violent",3),rep("Drug",3),rep("Low Level",3))
mods$year<-rep(c(1,2,3),4)
mods$Change<-(mods$Estimate-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100

write.csv(mods, file="results_ci.csv")


##-------------------------------------------------------------------------
## Main results
##-------------------------------------------------------------------------
m1<-glm.nb(all ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)
m2<-glm.nb(violent ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)
m3<-glm.nb(drug ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)
m4<-glm.nb(lowlevel ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)

mod1<-ses(m1)[52:54,]
mod2<-ses(m2)[52:54,]
mod3<-ses(m3)[52:54,]
mod4<-ses(m4)[52:54,]

mods<-rbind(mod1,mod2,mod3,mod4)
mods$Estimate<-exp(mods$Estimate)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$model<-c(rep("All",3),rep("Violent",3),rep("Drug",3),rep("Low Level",3))
mods$year<-rep(c(1,2,3),4)
mods$Change<-(mods$Estimate-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100

write.csv(mods, file="results_main.csv")

##-------------------------------------------------------------------------
## S3 Fig
##-------------------------------------------------------------------------

### SI Figs comparing Main, Coverage Indicator, PSM, PS, DOJ, Crime Adjusted, and State-level results
main <- data.frame(read.csv("results_main.csv", header=TRUE, stringsAsFactors =FALSE))
main$set<-"Main"
ci <- data.frame(read.csv("results_ci.csv", header=TRUE, stringsAsFactors =FALSE))
ci$set<-"Coverage Indicator"
psm <- data.frame(read.csv("results_psm.csv", header=TRUE, stringsAsFactors =FALSE))
psm$set<-"Propensity Score Matched"
state <- data.frame(read.csv("results_state.csv", header=TRUE, stringsAsFactors =FALSE))
state$set<-"State-Level"
crime <- data.frame(read.csv("results_crimeadj.csv", header=TRUE, stringsAsFactors =FALSE))
crime$set<-"Crime Adjusted"

plot<-rbind(main,ci,psm,state,crime)
plot$test<-plot$set

#png(file="S3fig_RR1.png",width=1200, height=700)
p1<-ggplot(plot[plot$model=="All" & plot$year==1,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("All")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "horizontal",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p2<-ggplot(plot[plot$model=="Violent" & plot$year==1,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Violent")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p3<-ggplot(plot[plot$model=="Drug" & plot$year==1,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Drug")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p4<-ggplot(plot[plot$model=="Low Level" & plot$year==1,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Low Level")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p<-plot_grid(p1+ theme(legend.position="none"),
             p2+ theme(legend.position="none"),
             p3+ theme(legend.position="none"), 
             p4+ theme(legend.position="none"), nrow=1,  align="vh")

legend_b <- get_legend(
  p1 +
    guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom", legend.text=element_text(size=18)))

p5<-plot_grid(p,legend_b, ncol=1,rel_heights = c(1,.1))
print(p5)
#dev.off()

#png(file="S3fig_RR2.png",width=1200, height=700)
p1<-ggplot(plot[plot$model=="All" & plot$year==2,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("All")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "horizontal",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p2<-ggplot(plot[plot$model=="Violent" & plot$year==2,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Violent")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p3<-ggplot(plot[plot$model=="Drug" & plot$year==2,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Drug")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p4<-ggplot(plot[plot$model=="Low Level" & plot$year==2,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Low Level")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p<-plot_grid(p1+ theme(legend.position="none"),
             p2+ theme(legend.position="none"),
             p3+ theme(legend.position="none"), 
             p4+ theme(legend.position="none"), nrow=1,  align="vh")

legend_b <- get_legend(
  p1 +
    guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom", legend.text=element_text(size=18)))

p5<-plot_grid(p,legend_b, ncol=1,rel_heights = c(1,.1))
print(p5)
#dev.off()

#png(file="S3fig_RR3.png",width=1200, height=700)
p1<-ggplot(plot[plot$model=="All" & plot$year==3,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("All")+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "horizontal",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p2<-ggplot(plot[plot$model=="Violent" & plot$year==3,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Violent")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p3<-ggplot(plot[plot$model=="Drug" & plot$year==3,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Drug")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p4<-ggplot(plot[plot$model=="Low Level" & plot$year==3,], aes(x=set, y=Change,col=test))+
  geom_point(size=3)+
  geom_line()+
  geom_errorbar(aes(ymin=CLL, ymax=CUL), width=0.2)+
  coord_cartesian(ylim=c(-51,30))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ggtitle("Low Level")+
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black", size=20),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), 
        axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(), 
        axis.title.y=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p<-plot_grid(p1+ theme(legend.position="none"),
             p2+ theme(legend.position="none"),
             p3+ theme(legend.position="none"), 
             p4+ theme(legend.position="none"), nrow=1,  align="vh")

legend_b <- get_legend(
  p1 +
    guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom", legend.text=element_text(size=18)))

p5<-plot_grid(p,legend_b, ncol=1,rel_heights = c(1,.1))
print(p5)
#dev.off()

##---------------------------------------------------------------------------
## Region of Common Support S5 Fig
##---------------------------------------------------------------------------

labs <- paste("ACA Expansion:", c("Expansion", "Non-Expansion"))

#png("S5fig.png")
prs_df %>%
  mutate(treatment = ifelse(treatment == 1, labs[1], labs[2])) %>%
  ggplot(aes(x = pr_score)) +
  geom_histogram(color = "white") +
  facet_wrap(~treatment) +
  xlab("Probability of Medicaid Expansion") +
  theme_bw()
#dev.off()

##---------------------------------------------------------------------------
## Coverage Indicator S6 Fig
##---------------------------------------------------------------------------

#pdf("S6fig.pdf")
p<-ggplot(analysis, aes(x=coverage_indicator)) + geom_histogram(color="black", fill="white")+ylab("Frequency")+xlab("Coverage Indicator")
print(p)
#dev.off()
